%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program: process_output.m
% By: Stephie Fried, Kevan, Novan, William Peterman
% Date: Spring 2024
% Purpose: Processes fortran output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
close all;
clear all;

bpath = 'C:\Users\l1sxf17\Dropbox (FRB SF)\Research\Will\equity_efficiency_project\Replication';
datapath = 'C:\Users\l1sxf17\Dropbox (FRB SF)\Research\Will\equity_efficiency_project\Replication\Fortran_output';
cd(bpath);
%% Baseline values
% Parameters 
theta1 = 2; 
J = 81;
n = 0.011; 

%Import population weights
pop_temp1=xlsread('other/calibration_conesa.xls','Sheet3');
for i = 1 : J-1
    survive(i)=pop_temp1(i+1)/pop_temp1(i); %survive(i) = prob survive to i+1 conditional on being alive at i
end
pop=zeros(J,2);
pop(1,1)=1;
for i = 2:J-1
    pop(i,1)=survive(i-1)*pop(i-1,1)/(1+n);
end
pop(:,2)=pop(:,1)./sum(pop(:,1)); %population share with growth
popYoung =pop(1:46,1)./sum(pop(1:46,1)); %pop share with growth of working
popOld =pop(47:81,1)./sum(pop(47:81,1)); %pop share with growth retirees
popweights  = pop(:,2); 

%Import data
M = importdata('Fortran_output/baseline/aggregates_bigger.dat');

w_base = M(1); r_base = M(2); K_base = M(3); N_base = M(4); lambda0_base = M(6); lambdak_base = M(12);
prob_base = M(13);  giniw_base = M(18); H_base = M(21); Ep_base = M(27); pe_base = M(16);
Ec_base = M(22); Ky_base = M(23); Ke_base = M(25); Ny_base = M(24); Ne_base = M(26);  E_base = M(11);
G_base = M(7); tauss_base = M(8); C_base = M(15); tr_base = M(5); ss = M(9); ctilde_base = M(30); lambda1_base = M(13);
avg_earnings_base = M(17); alpha1 = M(33); alpha2=M(34); theta = M(35); delta = M(36); A = M(37); Ae = M(38); 
gamma = M(39); ebar = M(40); beta = M(41); chi = M(42); theta2 = M(43); lambda1 = M(45); 
lambdak = M(46); r = M(2); ginii_base=M(32); 

I_base = delta*K_base;
R_base = r_base + delta;
after_tax_r_base = (1-lambdak_base)*r_base;
Y_base = Ky_base.^.3.*Ny_base.^(1-.3-.03).*Ep_base.^.03;
resource_base = Y_base - C_base  - G_base - (n+delta)*K_base;
taue = 0.482*pe_base;
gdp_base = Y_base + pe_base*Ec_base;

%indicies
w_ind = 1;
r_ind =  2;
K_ind = 3;
N_ind = 4;
E_ind = 11;
giniw_ind = 7; %welfare
cev_ind = 108; %CEV on non-energy consumption
Ky_ind = 80;
Ny_ind = 81;
Ep_ind =84;

%Tax parameters 
lambda1_ind = 14; %progressivity parameter
lump_ind = 16;
lambdak_ind = 12;
lambda0_ind = 13; %level parameter
 
%Avearge and marginal labor tax rates
fm =[0.1:.1:2]; %fraction of mean labor income
avg_rate = 1 - lambda0_base*fm.^(-lambda1_base);
yhvec =[0.1:0.01:10]*avg_earnings_base;
marg_rate = 1 - lambda0_base*(1-lambda1_base)*yhvec.^(-lambda1_base)*avg_earnings_base^lambda1_base; 

after_tax_yh1 = lambda0_base*yhvec.*(yhvec./avg_earnings_base).^(-lambda1_base); %baseline
after_tax_yh2 = 0.833091701*yhvec.*(yhvec./avg_earnings_base).^(-lambda1_base);  %Average labor-tax rebate

%% Read in fortran results
clc
cd(datapath);
SK = importdata('simple_capital/simple_capital.dat'); %capital tax rebate
SL = importdata('simple_labor/simple_labor.dat'); %labor tax rebate
SLS = importdata('simple_lump/simple_lump.dat'); %lump-sum tax rebate
OPT = importdata('optimal/aggregate_optimal.dat'); %optimal 
OCEAN = importdata('throw_away/throw_away.dat'); %ocean
PROG = importdata('prog/aggregates.dat'); %optimal progressivity exp 
MEID =  importdata('income_dependent/aggregates_income_maxprog_final.dat'); %Simple Income-dependent rebate, maximize equity
Heaven_mat = importdata('heaven/revenue_from_heaven_very_refined.dat'); %Need to choose optimal from heaven results
[temp, ind] = max(Heaven_mat(:, cev_ind));
HEAVEN = Heaven_mat(ind, :);
big_tax_mat =  importdata('tax_size/taue125_2.dat'); %Need to choose optimal
[temp, ind] = max(big_tax_mat(:, 18));
BIGTAX = big_tax_mat(ind, :);
small_tax_mat =  importdata('tax_size/taue075_2.dat'); %Need to choose optimal
[temp, ind] = max(small_tax_mat(:, 18));
SMALLTAX = small_tax_mat(ind, :);
ebarzero_mat =  importdata('ebar/ebar_zero_tax_refined3.dat'); %Need to choose optimal
[temp, ind] = max(ebarzero_mat(:, 18));
EBARZERO = ebarzero_mat(ind, :);
ebar2x_mat =  importdata('ebar/doubleebar_tax_refined1.dat'); %Need to choose optimal
[temp, ind] = max(ebar2x_mat(:, 18));
EBAR2X= ebar2x_mat(ind, :);
cd(bpath);

policies = [OPT', SK', SL', SLS', PROG', MEID'];  %Imposes the result that the welfarem maximizing income-dependent rebate is the lump-sum rebate
output = policies(Ky_ind, :).^alpha1.*policies(Ny_ind, :).^(1-alpha1 -0.03).*policies(Ep_ind).^0.03;

save('Programs/Matlab/aggregates');



